rm(list=ls())

# load packages
pacman::p_load(lubridate, tidyverse, haven, 
               gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra,readxl)

setwd('put_your_wd_here')

# shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date <- mdy(allshootings$date)
top10 <- allshootings$shooting[allshootings$top10==1]
filloutdates <- function(x){
  states <- unique(x$contributor_state)
  daterange <- seq(range(x$donationdate)[1],range(x$donationdate)[2], by='day')
  out <- data.frame(expand.grid(states, daterange))
  out <- out %>% 
    arrange(Var1) %>% 
    rename(contributor_state = Var1, donationdate = Var2)
  x <- left_join(
    out, x, by=c('donationdate','contributor_state')
  ) 
  x$donationcount[is.na(x$donationcount)] <- 0
  x$contribution_receipt_amount[is.na(x$contribution_receipt_amount)] <- 0
  x$Group <- x$Group[1]
  return(x)
}

nra <- read_stata('datasets/NRADailyStateDonations.dta')
nra <- filloutdates(nra)

filtDates <- function(df){
  # set date object
  df <- nra
  df$date <- ymd(df$donationdate)
  
  # set min and max dates for shootings
  mindate <- min(df$donationdate) + months(2)
  maxdate <- max(df$donationdate) - months(2)
  
  # subset to shootings that fit data datarange
  dates <- allshootings %>%
    filter(date >= mindate & date <= maxdate & include==1)
  return(dates)
}

nradates <- filtDates(nra)

makeEstRDD <- function(out, treat, name, org, time=2){
  
  # cut pre-post dates
  min <- treat - months(time) 
  max <- treat + months(time) 
  
  # set running day
  out$running_day <- as.numeric(out$donationdate - treat)
  
  # estimate sd
  sd_est <- sd(out$donationcount,na.rm=T)
  house_rdd <- rdd_data(y=out$donationcount, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}

states <- read_excel('datasets/shooting_states.xlsx')
nradates <- left_join(nradates, states[,c('shooting','state','state2','state3')], by='shooting')
nra_results_instate <- vector('list', nrow(nradates))
nra_results_outstate <- vector('list', nrow(nradates))

for(i in 1:length(nra_results_instate)){
  
  # subset to in-state data
  nradatinstate <- nra[nra$contributor_state %in% c(nradates$state[i], nradates$state2[i], nradates$state3[i]),]
  
  # group by date and sum up count
  nradatinstate <- nradatinstate %>%
    group_by(donationdate) %>%
    summarise(donationcount = sum(donationcount, na.rm=T))
  
  # out of state data
  nradatoutstate <- nra[!nra$contributor_state %in% c(nradates$state[i], nradates$state2[i], nradates$state3[i]),]
  nradatoutstate <- nradatoutstate %>%
    group_by(donationdate) %>%
    summarise(donationcount = sum(donationcount, na.rm=T))
  
  # results
  nra_results_instate[[i]] <- tryCatch({makeEstRDD(nradatinstate, nradates$date[i], 
                                       name = nradates$shooting[i], 
                                       org='NRA PAC')}, error=function(e){
                                         return(NA)})
  nra_results_outstate[[i]] <- makeEstRDD(nradatoutstate, nradates$date[i], 
                                          name = nradates$shooting[i], org='NRA PAC')
  print(i)
}

# collapse results
nra_results_instate <- nra_results_instate %>% discard(~ is.na(.x[1]))
nra_results_instate <- bind_rows(nra_results_instate) %>%
  mutate(results='In State')
nra_results_outstate <- bind_rows(nra_results_outstate) %>%
  mutate(results='Out of State')

datforplotnra_number <- expand.grid(
  shooting = unique(allshootings$shooting[allshootings$include==1]),
  results = c('In State','Out of State'),
  outcome = c('NRA PAC')
) %>%
  left_join(bind_rows(nra_results_instate, nra_results_outstate), by=c(
    'shooting','results','outcome'
  )) 


calcMeta <- function(data, name, outcome, shooting, group_var){
  m.gen <- meta::metagen(TE = Coef,
                         seTE = se,
                         studlab = shooting,
                         data = data,
                         sm = "SMD",
                         fixed = FALSE,
                         random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    group_var = group_var,
                    statsig=factor(1))
  return(results)
}

m1 <- calcMeta(nra_results_instate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (In State)',outcome='NRA',
         shooting='test',group_var = 'Number')
         
m2 <- calcMeta(nra_results_outstate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (Out State)',outcome='NRA',
         shooting='test',group_var = 'Number')


   

giffords <- read_stata('datasets/GiffordsPACDailyStateDonations.dta')
giffords <- filloutdates(giffords)
giffordsdates <- filtDates(giffords)

giffordsdates <- left_join(giffordsdates, states[,c('shooting','state','state2','state3')], by='shooting')
giff_results_instate <- vector('list', nrow(giffordsdates))
giff_results_outstate <- vector('list', nrow(giffordsdates))

for(i in 1:length(giff_results_instate)){
  giffinstatedat <- giffords[giffords$contributor_state %in% c(giffordsdates$state[i], giffordsdates$state2[i], giffordsdates$state3[i]),]
  giffinstatedat <- giffinstatedat %>%
    group_by(donationdate) %>%
    summarise(donationcount = sum(donationcount, na.rm=T))
  giffoutstate <- giffords[!giffords$contributor_state %in% c(giffordsdates$state[i], giffordsdates$state2[i], giffordsdates$state3[i]),]
  giffoutstate <- giffoutstate %>%
    group_by(donationdate) %>%
    summarise(donationcount = sum(donationcount, na.rm=T))
  giff_results_instate[[i]] <- tryCatch({makeEstRDD(giffinstatedat, giffordsdates$date[i], 
                                                   name = giffordsdates$shooting[i], 
                                                   org='Giffords PAC')}, error=function(e){
                                                     return(NA)})
  giff_results_outstate[[i]] <- makeEstRDD(giffoutstate, giffordsdates$date[i], 
                                          name = giffordsdates$shooting[i], org='Giffords PAC')
  print(i)
}

giff_results_instate <- giff_results_instate %>% discard(~ is.na(.x[1]))
giff_results_instate <- bind_rows(giff_results_instate) %>%
  mutate(results='In State')
giff_results_outstate <- bind_rows(giff_results_outstate) %>%
  mutate(results='Out of State')

datforplotgiff_number <- expand.grid(
  shooting = unique(allshootings$shooting[allshootings$include==1]),
  results = c('In State','Out of State'),
  outcome = c('Giffords PAC')
) %>%
  left_join(bind_rows(giff_results_instate, giff_results_outstate), by=c(
    'shooting','results','outcome'
  )) 

m3 <- calcMeta(giff_results_instate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (In State)',outcome='Giffords',
         shooting='test',group_var = 'Number')

m4 <- calcMeta(giff_results_outstate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (Out State)',outcome='Giffords',
         shooting='test',group_var = 'Number')

# AMOUNT

makeEstRDD <- function(out, treat, name, org, time=2){
  
  min <- treat - months(time) 
  max <- treat + months(time) 
  out$running_day <- as.numeric(out$donationdate - treat)
  
  sd_est <- sd(out$contribution_receipt_amount,na.rm=T)
  house_rdd <- rdd_data(y=out$contribution_receipt_amount, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}

states <- read_excel('datasets/shooting_states.xlsx')

# nradates <- left_join(nradates, states[,c('shooting','state','state2','state3')], by='shooting')
nra_results_instate <- vector('list', nrow(nradates))
nra_results_outstate <- vector('list', nrow(nradates))

for(i in 1:length(nra_results_instate)){
  
  # in state data
  nradatinstate <- nra[nra$contributor_state %in% c(nradates$state[i], nradates$state2[i], nradates$state3[i]),]
  nradatinstate <- nradatinstate %>%
    group_by(donationdate) %>%
    summarise(contribution_receipt_amount = sum(contribution_receipt_amount, na.rm=T))
  
  # out of state data
  nradatoutstate <- nra[!nra$contributor_state %in% c(nradates$state[i], nradates$state2[i], nradates$state3[i]),]
  nradatoutstate <- nradatoutstate %>%
    group_by(donationdate) %>%
    summarise(contribution_receipt_amount = sum(contribution_receipt_amount, na.rm=T))
  
  # in state results
  nra_results_instate[[i]] <- tryCatch({makeEstRDD(nradatinstate, nradates$date[i], 
                                                   name = nradates$shooting[i], 
                                                   org='NRA PAC')}, error=function(e){
                                                     return(NA)})
  
  # out of state results
  nra_results_outstate[[i]] <- makeEstRDD(nradatoutstate, nradates$date[i], 
                                          name = nradates$shooting[i], 
                                          org='NRA PAC')
  print(i)
}

nra_results_instate <- nra_results_instate %>% discard(~ is.na(.x[1]))
nra_results_instate <- bind_rows(nra_results_instate) %>%
  mutate(results='In State')
nra_results_outstate <- bind_rows(nra_results_outstate) %>%
  mutate(results='Out of State')

datforplotnra_amount <- expand.grid(
  shooting = unique(allshootings$shooting[allshootings$include==1]),
  results = c('In State','Out of State'),
  outcome = c('NRA PAC')
) %>%
  left_join(bind_rows(nra_results_instate, nra_results_outstate), by=c(
    'shooting','results','outcome'
  )) 

m5 <- calcMeta(nra_results_instate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (In State)',outcome='NRA',
         shooting='test',group_var = 'Amount')

m6 <- calcMeta(nra_results_outstate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (Out State)',outcome='NRA',
         shooting='test',group_var = 'Amount')

giff_results_instate <- vector('list', nrow(giffordsdates))
giff_results_outstate <- vector('list', nrow(giffordsdates))

for(i in 1:length(giff_results_instate)){
  
  # in state 
  giffinstatedat <- giffords[giffords$contributor_state %in% c(giffordsdates$state[i], giffordsdates$state2[i], giffordsdates$state3[i]),]
  giffinstatedat <- giffinstatedat %>%
    group_by(donationdate) %>%
    summarise(contribution_receipt_amount = sum(contribution_receipt_amount, na.rm=T))
  
  # out of states
  giffoutstate <- giffords[!giffords$contributor_state %in% c(giffordsdates$state[i], giffordsdates$state2[i], giffordsdates$state3[i]),]
  giffoutstate <- giffoutstate %>%
    group_by(donationdate) %>%
    summarise(contribution_receipt_amount = sum(contribution_receipt_amount, na.rm=T))
  
  # run models
  giff_results_instate[[i]] <- tryCatch({makeEstRDD(giffinstatedat, giffordsdates$date[i], 
                                                    name = giffordsdates$shooting[i], 
                                                    org='Giffords PAC')}, error=function(e){
                                                      return(NA)})
  giff_results_outstate[[i]] <- makeEstRDD(giffoutstate, giffordsdates$date[i], 
                                           name = giffordsdates$shooting[i], org='Giffords PAC')
  print(i)
}

giff_results_instate <- giff_results_instate %>% discard(~ is.na(.x[1]))
giff_results_instate <- bind_rows(giff_results_instate) %>%
  mutate(results='In State')
giff_results_outstate <- bind_rows(giff_results_outstate) %>%
  mutate(results='Out of State')

datforplotgiff_amount <- expand.grid(
  shooting = unique(allshootings$shooting[allshootings$include==1]),
  results = c('In State','Out of State'),
  outcome = c('Giffords PAC')
) %>%
  left_join(bind_rows(giff_results_instate, giff_results_outstate), by=c(
    'shooting','results','outcome'
  )) 

m7 <- calcMeta(giff_results_instate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (In State)',outcome='Giffords',
         shooting='test',group_var = 'Amount')

m8 <- calcMeta(giff_results_outstate %>% 
           mutate(se = (upper - Coef)/1.96),
         name='Meta Estimate (Out State)',outcome='Giffords',
         shooting='test',group_var = 'Amount')

datall <- bind_rows(datforplotnra_number %>% mutate(z = 'Number'),
          datforplotgiff_number %>% mutate(z = 'Number'),
          datforplotnra_amount %>% mutate(z = 'Amount'),
          datforplotgiff_amount %>% mutate(z = 'Amount'))
datall$shooting <- factor(datall$shooting, levels=rev(allshootings$shooting))

nra_plots <- datall %>%
  filter(outcome == 'NRA PAC') %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=results)) +
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_pointrange(position=position_dodge(width=0.5), fill='white') +
  coord_flip() + 
  facet_wrap(~z) +
  theme_bw() +
  theme(legend.position='bottom') +
  scale_shape_manual(values=c(21,23), name='') +
  labs(x='', y='')

giff_plot <- datall %>%
  filter(outcome == 'Giffords PAC') %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=results)) +
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_pointrange(position=position_dodge(width=0.5), fill='white') +
  coord_flip() + 
  facet_wrap(~z) +
  theme_bw() +
  theme(legend.position='bottom') +
  scale_shape_manual(values=c(21,23), name='') +
  labs(x='', y='')

nra_plots <- datall %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=results)) +
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_pointrange(position=position_dodge(width=0.5), fill='white') +
  coord_flip() + 
  facet_wrap(~outcome + z) +
  theme_bw() +
  theme(legend.position='bottom') +
  scale_shape_manual(values=c(21,23), name='') +
  labs(x='', y='')

meta_plot <- bind_rows(m1, m2, m3, m4, m5, m6, m7, m8) %>%
  mutate(shooting = gsub('Meta Estimate \\(|\\)','',shooting)) %>%
  ggplot(aes(x=group_var, y=Coef, ymin=lower, ymax=upper, shape=shooting)) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_pointrange(position=position_dodge(width=0.3), fill='white') +
  coord_flip() + 
  facet_wrap(~outcome) +
  theme_bw() +
  theme(legend.position='bottom') +
  scale_shape_manual(values=c(21,23), name='') +
  labs(x='', y='')
ggsave(meta_plot, width=5, height=3, 
       filename = 'figures/donations_state_RDITs_meta.pdf')


